R Markdown

chapter 3

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stringr)
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
## 
##     extract
library(spData)
world
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>
coffee_data
## # A tibble: 47 x 3
##    name_long                coffee_production_2016 coffee_production_2017
##    <chr>                                     <int>                  <int>
##  1 Angola                                       NA                     NA
##  2 Bolivia                                       3                      4
##  3 Brazil                                     3277                   2786
##  4 Burundi                                      37                     38
##  5 Cameroon                                      8                      6
##  6 Central African Republic                     NA                     NA
##  7 Congo, Dem. Rep. of                           4                     12
##  8 Colombia                                   1330                   1169
##  9 Costa Rica                                   28                     32
## 10 Côte d'Ivoire                               114                    130
## # … with 37 more rows
world_coffee = left_join(world, coffee_data) #左に右を結合(ここでは共通項name_longに依って結合.指定可能)
## Joining, by = "name_long"
class(world_coffee)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
world_coffee
## Simple feature collection with 177 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 13
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 4 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## #   coffee_production_2017 <int>
plot(world_coffee["coffee_production_2017"]) #地図情報とコーヒー情報を元にプロット

coffee_world = left_join(coffee_data, world)
## Joining, by = "name_long"
coffee_world
## # A tibble: 47 x 13
##    name_long     coffee_production… coffee_productio… iso_a2 continent region_un
##    <chr>                      <int>             <int> <chr>  <chr>     <chr>    
##  1 Angola                        NA                NA AO     Africa    Africa   
##  2 Bolivia                        3                 4 BO     South Am… Americas 
##  3 Brazil                      3277              2786 BR     South Am… Americas 
##  4 Burundi                       37                38 BI     Africa    Africa   
##  5 Cameroon                       8                 6 CM     Africa    Africa   
##  6 Central Afri…                 NA                NA CF     Africa    Africa   
##  7 Congo, Dem. …                  4                12 <NA>   <NA>      <NA>     
##  8 Colombia                    1330              1169 CO     South Am… Americas 
##  9 Costa Rica                    28                32 CR     North Am… Americas 
## 10 Côte d'Ivoire                114               130 CI     Africa    Africa   
## # … with 37 more rows, and 7 more variables: subregion <chr>, type <chr>,
## #   area_km2 <dbl>, pop <dbl>, lifeExp <dbl>, gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>
plot(coffee_world["coffee_production_2017"])

world_coffee_inner = inner_join(world, coffee_data) #left_joinしつつNA項を削る
## Joining, by = "name_long"
world_coffee_inner
## Simple feature collection with 45 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -117.1278 ymin: -33.76838 xmax: 156.02 ymax: 35.49401
## Geodetic CRS:  WGS 84
## # A tibble: 45 x 13
##    iso_a2 name_long  continent region_un subregion type  area_km2    pop lifeExp
##    <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>  <dbl>   <dbl>
##  1 TZ     Tanzania   Africa    Africa    Eastern … Sove…  932746. 5.22e7    64.2
##  2 PG     Papua New… Oceania   Oceania   Melanesia Sove…  464520. 7.76e6    65.2
##  3 ID     Indonesia  Asia      Asia      South-Ea… Sove… 1819251. 2.55e8    68.9
##  4 KE     Kenya      Africa    Africa    Eastern … Sove…  590837. 4.60e7    66.2
##  5 DO     Dominican… North Am… Americas  Caribbean Sove…   48158. 1.04e7    73.5
##  6 TL     Timor-Les… Asia      Asia      South-Ea… Sove…   14715. 1.21e6    68.3
##  7 MX     Mexico     North Am… Americas  Central … Sove… 1969480. 1.24e8    76.8
##  8 BR     Brazil     South Am… Americas  South Am… Sove… 8508557. 2.04e8    75.0
##  9 BO     Bolivia    South Am… Americas  South Am… Sove… 1085270. 1.06e7    68.4
## 10 PE     Peru       South Am… Americas  South Am… Sove… 1309700. 3.10e7    74.5
## # … with 35 more rows, and 4 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## #   coffee_production_2017 <int>
plot(world_coffee_inner["coffee_production_2017"])

world_new = world # do not overwrite our original data
world_new$pop_dens = world_new$pop / world_new$area_km2 #列の追加
world_new
## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 12
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##  * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, pop_dens <dbl>
world %>% 
  mutate(pop_dens = pop / area_km2) #格納してないのでworldデータ自体は変わっていない
## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 12
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##  * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, pop_dens <dbl>
world %>% transmute(pop_dens = pop / area_km2) %>% plot() #pop_densのみにする(sfクラスなのでgeomは残る)

elev = raster(nrows = 6, ncols = 6, res = 0.5,
              xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
              vals = 1:36) #左から右、その後下の行
plot(elev)

grain_order = c("clay", "silt", "sand")
grain_char = sample(grain_order, 36, replace = TRUE) #36個にランダムに3つをあてがう
grain_fact = factor(grain_char, levels = grain_order) #rasterが扱えるようにしてる
grain = raster(nrows = 6, ncols = 6, res = 0.5, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = grain_fact)
plot(grain)

elev[1,1]
## [1] 1
elev[1]
## [1] 1
r_stack = stack(elev, grain) #データのスタック
names(r_stack) = c("elev", "grain")
# three ways to extract a layer of a stack
raster::subset(r_stack, "elev")
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)
r_stack[["elev"]]
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)
r_stack$elev
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)
elev[]
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36
elev[1, 1:2] = 0

cellStats(elev, sd)
## [1] 10.67808
hist(elev)

chapter 4

library(sf)
library(raster)
library(dplyr)
library(spData)
nz
## Simple feature collection with 16 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1090144 ymin: 4748537 xmax: 2089533 ymax: 6191874
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## First 10 features:
##                 Name Island Land_area Population Median_income Sex_ratio
## 1          Northland  North 12500.561     175500         23400 0.9424532
## 2           Auckland  North  4941.573    1657200         29600 0.9442858
## 3            Waikato  North 23900.036     460100         27900 0.9520500
## 4      Bay of Plenty  North 12071.145     299900         26200 0.9280391
## 5           Gisborne  North  8385.827      48500         24400 0.9349734
## 6        Hawke's Bay  North 14137.524     164000         26100 0.9238375
## 7           Taranaki  North  7254.480     118000         29100 0.9569363
## 8  Manawatu-Wanganui  North 22220.608     234500         25000 0.9387734
## 9         Wellington  North  8048.553     513900         32700 0.9335524
## 10        West Coast  South 23245.456      32400         26900 1.0139072
##                              geom
## 1  MULTIPOLYGON (((1745493 600...
## 2  MULTIPOLYGON (((1803822 590...
## 3  MULTIPOLYGON (((1860345 585...
## 4  MULTIPOLYGON (((2049387 583...
## 5  MULTIPOLYGON (((2024489 567...
## 6  MULTIPOLYGON (((2024489 567...
## 7  MULTIPOLYGON (((1740438 571...
## 8  MULTIPOLYGON (((1866732 566...
## 9  MULTIPOLYGON (((1881590 548...
## 10 MULTIPOLYGON (((1557042 531...
canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
plot(canterbury_height)

library(tmap)
p_hpnz1 = tm_shape(nz) + tm_polygons(col = "white") +
  tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 0.25) +
  tm_layout(main.title = "High points in New Zealand", main.title.size = 1,
            bg.color = "lightblue")
p_hpnz2 = tm_shape(nz) + tm_polygons(col = "white") +
  tm_shape(canterbury) + tm_fill(col = "gray") + 
  tm_shape(canterbury_height) + tm_symbols(shape = 2, col = "red", size = 0.25) +
  tm_layout(main.title = "High points in Canterbury", main.title.size = 1,
            bg.color = "lightblue")
tmap_arrange(p_hpnz1, p_hpnz2, ncol = 2)

sel_sgbp = st_intersects(x = nz_height, y = canterbury)
class(sel_sgbp)
## [1] "sgbp" "list"
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]

a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a = st_sfc(a_poly)
# create a line
l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2))
l = st_sfc(l_line)
# create points
p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi = st_multipoint(x = p_matrix)
p = st_cast(st_sfc(p_multi), "POINT")

par(pty = "s")
plot(a, border = "red", col = "gray", axes = TRUE)
plot(l, add = TRUE)
plot(p, add = TRUE, lab = 1:4)
text(p_matrix[, 1] + 0.04, p_matrix[, 2] - 0.06, 1:4, cex = 1.3)

sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix
lengths(sel) > 0
## [1]  TRUE  TRUE FALSE  TRUE
set.seed(2018) # set seed for reproducibility
(bb_world = st_bbox(world)) # the world's bounds(世界の範囲)
##       xmin       ymin       xmax       ymax 
## -180.00000  -90.00000  180.00000   83.64513
random_df = tibble(
  x = runif(n = 10, min = bb_world[1], max = bb_world[3]),
  y = runif(n = 10, min = bb_world[2], max = bb_world[4])
) #世界の範囲内でランダム点10個
random_points = random_df %>% 
  st_as_sf(coords = c("x", "y")) %>% # set coordinates
  st_set_crs(4326) # set geographic CRS

world_random = world[random_points, ] #random_pointsを国土内(geom)に含む国が出てくる(4/10)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
nrow(world_random)
## [1] 4
random_joined = st_join(random_points, world["name_long"])
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
plot(random_joined)

plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")

cycle_hire_P = st_transform(cycle_hire, 27700) #座標系を27700系に移す
cycle_hire_osm_P = st_transform(cycle_hire_osm, 27700)
sel = st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20) #hireとhire_osmで距離20m以内ならselに入れる
summary(lengths(sel) > 0)
##    Mode   FALSE    TRUE 
## logical     304     438
z = st_join(cycle_hire_P, cycle_hire_osm_P,
            join = st_is_within_distance, dist = 20)
nrow(cycle_hire)
## [1] 742
plot(cycle_hire)

nrow(z)
## [1] 762
z
## Simple feature collection with 762 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 522502 ymin: 174408 xmax: 538733.2 ymax: 184421
## Projected CRS: OSGB 1936 / British National Grid
## First 10 features:
##    id             name.x             area nbikes nempty    osm_id
## 1   1       River Street      Clerkenwell      4     14 869697014
## 2   2 Phillimore Gardens       Kensington      2     34 885331201
## 3   3 Christopher Street Liverpool Street      0     32 920087626
## 4   4  St. Chad's Street     King's Cross      4     19 781506147
## 5   5     Sedding Street    Sloane Square     15     12      <NA>
## 6   6 Broadcasting House       Marylebone      0     18 839403852
## 7   7   Charlbert Street  St. John's Wood     15      0 839403885
## 8   8         Lodge Road  St. John's Wood      5     13 839403863
## 9   9     New Globe Walk         Bankside      3     16 835264541
## 10 10        Park Street         Bankside      1     17 848525385
##                            name.y capacity cyclestreets_id description
## 1                    River Street        9            <NA>        <NA>
## 2  Kensington, Phillimore Gardens       27            <NA>        <NA>
## 3              Christopher Street       NA            <NA>        <NA>
## 4                            <NA>       NA            <NA>        <NA>
## 5                            <NA>       NA            <NA>        <NA>
## 6              Broadcasting House        8            <NA>        <NA>
## 7                Charlbert Street        2            <NA>        <NA>
## 8                      Lodge Road        8            <NA>        <NA>
## 9        New Globe Walk, Bankside        9            <NA>        <NA>
## 10                    Park Street        8            <NA>        <NA>
##                     geometry
## 1  POINT (531203.5 182832.1)
## 2  POINT (525208.1 179391.9)
## 3  POINT (532985.8 182001.6)
## 4    POINT (530437.8 182912)
## 5      POINT (528051 178742)
## 6  POINT (528858.4 181542.9)
## 7    POINT (527159 183300.8)
## 8  POINT (527032.7 182634.6)
## 9    POINT (532205 180434.6)
## 10 POINT (532464.9 180284.3)
plot(nz)

plot(nz_height)

nz_avheight = aggregate(x = nz_height, by = nz, FUN = mean) #FUNの引数を変えることで様々な関数になる

library(tmap)
tm_shape(nz_avheight) +
  tm_fill("elevation", breaks = seq(27, 30, by = 0.5) * 1e2) +
  tm_borders()

plot(elev)

id = cellFromXY(elev, xy = c(0.1, 0.1)) #xyが(0.1,0.1)のインデックスをとる
elev[id]
## [1] 16
# the same as
raster::extract(elev, data.frame(x = 0.1, y = 0.1))
## [1] 16
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
              res = 0.3, vals = rep(1, 9))
plot(clip)

elev[clip]
## [1] 18 24
# create raster mask
rmask = elev 
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)
plot(rmask)

# spatial subsetting
elev[rmask, drop = FALSE]           # with [ operator
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 36  (min, max)
mask(elev, rmask)                   # with mask()
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 36  (min, max)
overlay(elev, rmask, fun = "max")   # with overlay
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 36  (min, max)
elev + elev #セルサイズが同じなら演算が可能
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 72  (min, max)
elev^2
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1296  (min, max)
log(elev)
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -Inf, 3.583519  (min, max)
elev > 5
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

chapter 5

library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
## 
## Attaching package: 'spDataLarge'
## The following object is masked _by_ '.GlobalEnv':
## 
##     random_points
us_states
## Simple feature collection with 49 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
## Geodetic CRS:  NAD83
## First 10 features:
##    GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
## 1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
## 2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
## 3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
## 4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
## 5     12     Florida    South 151052.01 [km^2]     18511620     19645772
## 6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
## 7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
## 8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
## 9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
## 10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
##                          geometry
## 1  MULTIPOLYGON (((-88.20006 3...
## 2  MULTIPOLYGON (((-114.7196 3...
## 3  MULTIPOLYGON (((-109.0501 4...
## 4  MULTIPOLYGON (((-73.48731 4...
## 5  MULTIPOLYGON (((-81.81169 2...
## 6  MULTIPOLYGON (((-85.60516 3...
## 7  MULTIPOLYGON (((-116.916 45...
## 8  MULTIPOLYGON (((-87.52404 4...
## 9  MULTIPOLYGON (((-102.0517 4...
## 10 MULTIPOLYGON (((-92.01783 2...
plot(us_states)

regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
plot(regions)

regions2 = us_states %>% group_by(REGION) %>%
  summarize(pop = sum(total_pop_15, na.rm = TRUE))
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west) #unionは中身を消してgeomのみにする?
## although coordinates are longitude/latitude, st_union assumes that they are planar
plot(us_west)

plot(us_west_union)

data("dem", package = "spDataLarge")
plot(dem)

dem
## class      : RasterLayer 
## dimensions : 117, 117, 13689  (nrow, ncol, ncell)
## resolution : 30.85, 30.85  (x, y)
## extent     : 794599.1, 798208.6, 8931775, 8935384  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=17 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : memory
## names      : dem 
## values     : 238, 1094  (min, max)
dem_agg = aggregate(dem, fact = 5, fun = mean) #解像度を下げる(1/5)
plot(dem_agg)

dem_agg
## class      : RasterLayer 
## dimensions : 24, 24, 576  (nrow, ncol, ncell)
## resolution : 154.25, 154.25  (x, y)
## extent     : 794599.1, 798301.1, 8931682, 8935384  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=17 +south +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : memory
## names      : dem 
## values     : 239.9, 1077.6  (min, max)
dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear") #解像度を上げる(5倍)
plot(dem_disagg)

identical(dem, dem_disagg) #完全には戻らない
## [1] FALSE
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge")) #標高データ
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge")) #zionという場所(srtmが示している場所)のデータ
## Reading layer `zion' from data source `/usr/local/lib/R/site-library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## Projected CRS: UTM Zone 12, Northern Hemisphere
zion = st_transform(zion, projection(srtm)) #2つのデータで座標系を統一
plot(srtm)

plot(zion)
## Warning: plotting the first 10 out of 11 attributes; use max.plot = 11 to plot
## all

srtm_cropped = crop(srtm, zion) #zionの地形全体を含む長方形でsrtmをクロップしている
plot(srtm_cropped)

srtm_masked = mask(srtm, zion) #zionの形に沿ってクロップ(マスク)
plot(srtm_masked)

data("zion_points", package = "spDataLarge")
zion_points$elevation = raster::extract(srtm, zion_points)
#raster::extract(srtm, zion_points, buffer = 1000) #バッファ込み.めちゃ重い

zion_nlcd = raster::extract(nlcd, zion, df = TRUE, factors = TRUE) 
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
## Warning in spTransform(xSP, CRSobj, ...): NULL target CRS comment, falling back
## to PROJ string
dplyr::select(zion_nlcd, ID, levels) %>% 
  tidyr::gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>% 
  tidyr::spread(value, n, fill = 0)
## # A tibble: 1 x 9
## # Groups:   ID, key [1]
##      ID key    Barren Cultivated Developed Forest Herbaceous Shrubland Wetlands
##   <dbl> <chr>   <dbl>      <dbl>     <dbl>  <dbl>      <dbl>     <dbl>    <dbl>
## 1     1 levels  98285         62      4205 298299        235    203701      679
cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4
## definition
ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
plot(ch_raster1)

ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = 1, fun = "count")
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
plot(ch_raster2)

ch_raster3 = rasterize(cycle_hire_osm_projected, raster_template, 
                       field = "capacity", fun = sum)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
plot(ch_raster3)